
//load IPCC scenario data 
use ../Intermediate/ipcc_scenarios_data, clear
drop Region Model

gen scen = "_current" if Scenario == "NGFS2_Current Policies"
replace scen = "_below2c" if Scenario == "CO_Bridge" 
drop Scenario 
drop  primaryenergy_coal_EJyr
rename secondaryenergy_elec_coal_EJyr coal_elec_EJ 
rename kyoto_gases_Mt_CO2eyr total_emissions_Mt_CO2eyr

//convert coal electricity to CO2e Emssions
gen coal_emissions_Mt_CO2eyr = coal_elec_EJ*0.228*1000
keep scen year coal_emissions_Mt_CO2eyr coal_power_capacity_GW
reshape wide coal_*, i(year) j(scen) string
tempfile ipcc 
save `ipcc'

clear
set obs 91
gen year = 2010+_n-1
tempfile time 
save `time'


//get hazard model sample to understand the overal plant level treatment effect
use ../Intermediate/plant_year_panel_clean, clear
stset year, failure(retired) id(plant_parent_id) origin(EstYr) 
stcox c.post2015##c.ban_intensity_max_sd, strata(country_i) vce(cl ParentID_GCEL)
gen insamp = e(sample) == 1
egen plant_in_cox_samp = max(insamp), by(TrackerID)
keep TrackerID ban_intensity_max_sd plant_in_cox_samp 
collapse (mean) ban_intensity_max_sd, by(TrackerID plant_in_cox_samp)
duplicates drop
tempfile coxsamp 
save `coxsamp'

use "../Intermediate/gcpt_clean.dta", clear

unique TrackerID // unique plant identifier 

gen ever_operated = ~missing(first_date_Operating)
replace ever_operated = 0 if missing(EstYr) & inlist(most_recent_status_2022,"Cancelled","Construction")
replace ever_operated = 0 if ~inlist(most_recent_status_2022,"Operating","Retired","Mothballed")


gen EstYr_imputed = EstYr if EstYr <= 2022
egen EstYr_location = mean(EstYr_imputed), by(TrackerLOC)
replace EstYr_location = floor(EstYr_location)
replace EstYr_imputed = EstYr_location if missing(EstYr) & ever_operated == 1
mdesc EstYr_imputed if ever_operated == 1
egen EstYr_province = mean(EstYr_imputed), by(PlantProvince)
replace EstYr_province = floor(EstYr_province)
replace EstYr_imputed = EstYr_province if missing(EstYr_imputed) & ever_operated == 1
mdesc EstYr_imputed if ever_operated == 1
replace EstYr_imputed = EstYr if missing(EstYr_imputed) & ever_operated == 0
gen future_plants = EstYr_imputed >= 2023 & ~missing(EstYr_imputed) & ~inlist(most_recent_status_2022,"Cancelled","Shelved","Retired")


//GENERATE SAMPLE OF PLANTS ACTIVE AS OF OR AFTER 2010
keep if (ever_operated == 1 & EstYr_imputed <= 2022 & inlist(most_recent_status_2022,"Operating","Retired","Mothballed")) | future_plants == 1


//survival analysis
gen dox = RetiredYr
replace dox = 2022 if missing(RetiredYr)
gen retired = ~missing(RetiredYr)
gen entered = max(EstYr_imputed,2014) 
stset dox, failure(retired) enter(time entered) origin(time EstYr_imputed)
encode PlantCountry, gen(country_i)

preserve
drop if future_plants == 1
sts graph, name(KMsurv, replace)
sts generate surv_km = s
streg, distribution(weibull)
keep surv_km _t 
duplicates drop 
tempfile km 
save `km'
clear
set obs 81
gen t0 = 0
gen t = _n 
gen st = 1
gen d = 1
stset t, failure(d) 
predict surv_weibull, surv
predict haz_weibull, haz
merge 1:1 _t using `km', keep(1 3)
line surv_weibull surv_km _t, scheme(white_tableau) ///
							  xtitle("Plant Age") ///
						      legend(pos(6) cols(2) order(2 "Kaplan-Meier Survival Estimate" 1 "Weibull Model Fit"))


//hazard rate 
keep surv_weibull haz_weibull _t 
rename _t t
tempfile weibull 
save `weibull'
restore

merge 1:1 TrackerID using `coxsamp', keep(1 3) nogen
gen ban_intensity = ban_intensity_max_sd
replace ban_intensity = 0 if missing(ban_intensity_max_sd)

keep TrackerID EstYr_imputed future_plants RetiredYr CapacityMW Subregion Coaltype HeatrateBtuperkWh Capacityfactor AnnualCO2milliontonnesann Remainingplantlifetimeyears PlantCountry
rename EstYr_imputed EstYr 
order TrackerID EstYr RetiredYr

local haz_mult_nopol  = 1/exp(0.366*0.497)  //multiply baseline hazard by this to get no policy hazard 
local haz_mult_bigpol = exp(0.366*1.729)/exp(0.366*0.497)  //multiply baseline hazard by this to get bigger policy hazard


//make panel from 2010
cross using `time'
order year 

egen plantid = group(TrackerID)
tsset plantid year
gen retired = year>=RetiredYr
gen entered = max(2014,EstYr)
gen exit = min(RetiredYr,2022)
stset year, failure(retired) id(plantid) origin(time EstYr) enter(time entered) exit(time exit)
gen t = year - EstYr  //negative years means plant was not yet established
merge m:1 t using `weibull'


tsset plantid year

//actual data
gen data_operating = year < RetiredYr & year >= EstYr
gen data_retired = year >= RetiredYr 

//use data for operating status before 2015 or through plant establishment
gen baseline_p_alive = data_operating if year <= max(EstYr,2015)
replace baseline_p_alive = 0 if year > 2015 & RetiredYr <= 2015 //already retired plants stay retired
//after 2015 or estyear, start applying hazard
bys plantid (year): replace baseline_p_alive = baseline_p_alive[_n-1]*(1-haz_weibull) if year > max(EstYr,2015)

gen nopolicy_p_alive = data_operating if year <= max(EstYr,2015)
replace nopolicy_p_alive = 0 if year > 2015 & RetiredYr <= 2015 //already retired plants stay retired
//after 2015 or estyear, start applying hazard
gen haz_nopol = haz_weibull*`haz_mult_nopol'
bys plantid (year): replace nopolicy_p_alive = nopolicy_p_alive[_n-1]*(1-haz_nopol) if year > max(EstYr,2015)


gen bigpol_p_alive = data_operating if year <= max(EstYr,2015)
replace bigpol_p_alive = 0 if year > 2015 & RetiredYr <= 2015 //already retired plants stay retired
//after 2015 or estyear, start applying hazard
gen haz_bigpol = haz_weibull*`haz_mult_bigpol'
bys plantid (year): replace bigpol_p_alive = bigpol_p_alive[_n-1]*(1-haz_bigpol) if year > max(EstYr,2015)

//make a COAL emissions over time plot 
gen data_co2eGT = data_operating*AnnualCO2milliontonnesann/1e3
gen baseline_co2eGT = baseline_p_alive*AnnualCO2milliontonnesann/1e3
gen nopolicy_co2eGT = nopolicy_p_alive*AnnualCO2milliontonnesann/1e3
gen bigpolicy_co2eGT = bigpol_p_alive*AnnualCO2milliontonnesann/1e3

gen data_capacityGW = data_operating*CapacityMW/1e3
gen baseline_capacityGW = baseline_p_alive*CapacityMW/1e3
gen nopolicy_capacityGW = nopolicy_p_alive*CapacityMW/1e3

//make GE EMISSIONS PLOT 
tempfile pre_collapse 
save `pre_collapse'

use `pre_collapse', clear
collapse (sum) CapacityMW data_co2eGT baseline_co2eGT nopolicy_co2eGT bigpolicy_co2eGT *_capacityGW, by(year)
label var nopolicy_co2eGT "No Policies"
label var baseline_co2eGT "Current Policies"
label var bigpolicy_co2eGT "Expanded Policies"

label var nopolicy_capacityGW "No Policies"
label var baseline_capacityGW "Current Policies"

label var year ""

merge 1:1 year using `ipcc'
drop if year < 2010
drop _merge
gen zero = 0

gen below2c_coal_co2eGT = coal_emissions_Mt_CO2eyr_below2c/1000
gen bau_coal_co2eGT = coal_emissions_Mt_CO2eyr_current/1000

label var below2c_coal_co2eGT "Below 2C Target"
label var bau_coal_co2eGT "IPCC Current Path"

save "../Intermediate/Figure_8", replace
 


